## Load packages and set wd ----

if (!require(QCA)) {install.packages("QCA"); require(QCA)} # Generates truth-table-style tables
if (!require(egg)) {install.packages("egg"); require(egg)} # Puts multiple ggplots into one panel
if (!require(lme4)) {install.packages("lme4"); require(lme4)}
if (!require(pastecs)) {install.packages("pastecs"); require(pastecs)}
if (!require(car)) {install.packages("car"); require(car)}
if (!require(grDevices)) {install.packages("grDevices"); require(grDevices)}
if (!require(effects)) {install.packages("effects"); require(effects)}
if (!require(optimx)) {install.packages("optimx"); require(optimx)}
if (!require(emmeans)) {install.packages("emmeans"); require(emmeans)}
if (!require(stringr)) {install.packages("stringr"); require(stringr)}
if (!require(rstudioapi)) {install.packages("rstudioapi"); require(rstudioapi)}
if (!require(WRS2)) {install.packages("WRS2"); require(WRS2)}
if (!require(broom)) {install.packages("broom"); require(broom)}
if (!require(tidyverse)) {install.packages("tidyverse"); require(tidyverse)}


current_path <- getActiveDocumentContext()$path    ## gets path for current R script
setwd(dirname(current_path))                       ## sets dir to R script path

cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")

## Load data for analyzing children's responses to modal questions ----

mq.df.long <- read.csv("mq.df.long.csv")

n.by.age <- mq.df.long %>% # n per age group; this will be used for plotting later
  group_by(age.group) %>%
  summarize(count = n_distinct(id)) %>%
  mutate(count.exp = paste("n = ", count, sep = ""))

# Let's have a look at the responses kids gave
length(unique(mq.df.long$response)) # 135 different response types (to 2304 questions)
table(mq.df.long$response)

# Let's code the responses into categories

## Principles for categorization:
# When a participant answers + elaborates, the answer determines the code
# When a participant answers + corrects, the correction determines the code
# Necessity modals were categorized as "yes"(They were rare).
# Negated necessity modals were categorized as "possibility modals" (also rare).
# When it was unclear whether a response was an elaboration or a correction, we used context to evaluate whether we
# thought the participant understood the question, and if so, coded their response in a way that reflected comprehension.

# Sometimes participants seemed to answer and then make a guess; we coded their answer.

mq.df.long$response.recode <- mq.df.long$response
mq.df.long$response.recode[mq.df.long$response.recode == "yes " |
                             mq.df.long$response.recode == "is.this.a.yes.no.question...yes" |
                             mq.df.long$response.recode == "i.think.so" |
                             mq.df.long$response.recode == "i.think.so." |
                             mq.df.long$response.recode == "im.pretty.sure.it.has.to" | # Necessity modal
                             mq.df.long$response.recode == "it.does" |
                             mq.df.long$response.recode == "maybe.yes" | # Counting this as a correction
                             mq.df.long$response.recode == "mm.yes" |
                             mq.df.long$response.recode == "no..yes" |
                             mq.df.long$response.recode == "no.i.think.so" |
                             mq.df.long$response.recode == "no.uh.yes" |
                             mq.df.long$response.recode == "no.yes" |
                             mq.df.long$response.recode == "sure" | 
                             mq.df.long$response.recode == "sure.if.it.wants.to" | # Necessity modal + desire verb
                             mq.df.long$response.recode == "uh.huh" |
                             mq.df.long$response.recode == "um.yes" |
                             mq.df.long$response.recode == "yeah" |
                             mq.df.long$response.recode == "yes.because.again.it.has.a.choice" | # Response + elaboration
                             mq.df.long$response.recode == "yes.because.its.the.only.way.the.marble.can.go" |
                             mq.df.long$response.recode == "yes.but.it.might.come.out.the.other.side.too" |
                             mq.df.long$response.recode == "yes.but.theres.the.other.side.too" |
                             mq.df.long$response.recode == "yes.it.can" | 
                             mq.df.long$response.recode == "yes.it.can.because.its.connected.so.it.has.a.choice.to.come.out.there" |
                             mq.df.long$response.recode == "yes.it.can.come.out.either.one" |
                             mq.df.long$response.recode == "yes.it.could.be.either.one" |
                             mq.df.long$response.recode == "yes.it.has.to" | # Yes + necessity modal
                             mq.df.long$response.recode == "yes.it.wil.because.thats.the.only.way.its.allowed" |
                             mq.df.long$response.recode == "yes.it.will" | # Yes + necessity modal
                             mq.df.long$response.recode == "yes.its.not.a.marble.thats.magic" |
                             mq.df.long$response.recode == "yes.maybe" |
                             mq.df.long$response.recode == "yes.sometimes" |
                             mq.df.long$response.recode == "yes.that's.the.only.way"] <- "yes"

mq.df.long$response.recode[mq.df.long$response.recode == "actually.no" |
                             mq.df.long$response.recode == "does.it.have.to.come.out.there.no" |
                             mq.df.long$response.recode == "does.not.go.there" |
                             mq.df.long$response.recode == "i.dont.think.so" |
                             mq.df.long$response.recode == "i.dont.think.so.either" |
                             mq.df.long$response.recode == "i.think.no" |
                             mq.df.long$response.recode == "it.wont.be.able.to.no" |
                             mq.df.long$response.recode == "maybe.no" | # Possibility modal plus correction
                             mq.df.long$response.recode == "no.because.it.cant" |
                             mq.df.long$response.recode == "no.because.its.not.connected" |
                             mq.df.long$response.recode == "no.because.its.not.connected.if.they.were.it.might" |
                             mq.df.long$response.recode == "no.because.that's.different" |
                             mq.df.long$response.recode == "no.because.they.don't.attach" |
                             mq.df.long$response.recode == "no.because.theyre.not.connected" |
                             mq.df.long$response.recode == "no.but.it.could" | 
                             mq.df.long$response.recode == "no.cuz.again.theyre.not.connected" |
                             mq.df.long$response.recode == "no.i.picked.the.other.one" |
                             mq.df.long$response.recode == "no.impossible" | # Elaboration?
                             mq.df.long$response.recode == "no.impossible.it.could" | # Elaboration?
                             mq.df.long$response.recode == "no.it.can.come.out.the.other.way" |
                             mq.df.long$response.recode == "no.it.can't" | # Elaboration?
                             mq.df.long$response.recode == "no.it.cant" | # Elaboration?
                             mq.df.long$response.recode == "no.it.cannot" |
                             mq.df.long$response.recode == "no.it.does.not" |
                             mq.df.long$response.recode == "no.it.doesnt.have.to.come.out.there" |
                             mq.df.long$response.recode == "no.it.wont" |
                             mq.df.long$response.recode == "no.it.will.not" |
                             mq.df.long$response.recode == "no.it's.not.above.it." |
                             mq.df.long$response.recode == "no.it's.not.attached" |
                             mq.df.long$response.recode == "no.not.necessarily" |
                             mq.df.long$response.recode == "no.thats.not.an.option" |
                             mq.df.long$response.recode == "nope" |
                             mq.df.long$response.recode == "no.that.way.is.different." |
                             mq.df.long$response.recode == "no.they're.not.connected"]<- "no"

mq.df.long$response.recode[mq.df.long$response.recode == "cannot" |
                             mq.df.long$response.recode == "impossible" |
                             mq.df.long$response.recode == "it.can't" |
                             mq.df.long$response.recode == "it.can't.at.all" |
                             mq.df.long$response.recode == "it.cannot.at.all" |
                             mq.df.long$response.recode == "it.cannot" |
                             mq.df.long$response.recode == "it.cant" |
                             mq.df.long$response.recode == "it.cant.no"] <- "negated.poss.modal" 

mq.df.long$response.recode[mq.df.long$response.recode == "can.i.say.maybe" |
                             mq.df.long$response.recode == "could.fifty.fifty.chance" |
                             mq.df.long$response.recode == "doesn't.hafta" | # Negated necessity modal
                             mq.df.long$response.recode == "i.does.not.hafta.cause.that.one.is.different" | # Negated necessity modal
                             mq.df.long$response.recode == "i.dont.know.it.could" |
                             mq.df.long$response.recode == "i.dont.know.it.might" |
                             mq.df.long$response.recode == "im.not.sure.it.might" |
                             mq.df.long$response.recode == "is.maybe.an.option" |
                             mq.df.long$response.recode == "it.can.but.it.might.not...no" | 
                             mq.df.long$response.recode == "it.could" |
                             mq.df.long$response.recode == "it.could " |
                             mq.df.long$response.recode == "it.could.i.dont.know" | 
                             mq.df.long$response.recode == "it.can" |
                             mq.df.long$response.recode == "it.can.because.that.one.has.two.slide" |
                             mq.df.long$response.recode == "it.can.come.out.there" |
                             mq.df.long$response.recode == "it.can.come.out.two" |
                             mq.df.long$response.recode == "it.can.it.does.not.have.to" | # Possibility modal + elaboration with negated necessity modal
                             mq.df.long$response.recode == "it.could " |
                             mq.df.long$response.recode == "it.could" |
                             mq.df.long$response.recode == "it.could.be.either.one" |
                             mq.df.long$response.recode == "it.could.come.out.the.other.one.on.that.one.too" |
                             mq.df.long$response.recode == "it.could.come.out.there.or.there" |
                             mq.df.long$response.recode == "it.could.go.out.the.other.one.you.pointed.to" |
                             mq.df.long$response.recode == "it.could.if.it.rolled.but.it.doesnt.have.to" | 
                             mq.df.long$response.recode == "it.does.not.have.to.no" | # Negated necessity modal (elaborated)
                             mq.df.long$response.recode == "it.doesnt.have.to.but.it.can" | # Negated necessity modal (elaborated)
                             mq.df.long$response.recode == "it.doesnt.have.to.but.it.could" | # Negated necessity modal (elaborated)
                             mq.df.long$response.recode == "it.might" |
                             mq.df.long$response.recode == "it.might.i.think.one.of.those.two.ones.it.might.come.out.each" |
                             mq.df.long$response.recode == "it's.one.of.it's.options" | 
                             mq.df.long$response.recode == "maybe" |
                             mq.df.long$response.recode == "maybe.we.don't.know" | 
                             mq.df.long$response.recode == "maybe.so.yes" | # The 'so' makes this seem like an elaboration, not a correction
                             mq.df.long$response.recode == "maybe.yes.i.think.ill.settle.on.yes" | # Adult; definitely sees the possibilities but also interprets the q as a prompt to guess
                             mq.df.long$response.recode == "might.fifty.fifty.chance" |
                             mq.df.long$response.recode == "no.uh.yeah.maybe" | # Correction; possibility modal
                             mq.df.long$response.recode == "or.maybe.the.other.side" |
                             mq.df.long$response.recode == "or.there" | # We read this as "It will come out there or there", and analyze disjunction as a possibility modal
                             mq.df.long$response.recode == "possibly.yes" | # Video checked; coded as elaboration on the basis of prosody
                             mq.df.long$response.recode == "same.thing.like.maybe.yes" | # Adult; definitely sees the possibilities but interprets the q as a prompt to guess
                             mq.df.long$response.recode == "same.thing.theres.no.certainty.so.it.could.but.i.guess.no" | # Adult; uses a possibility modal; understood modal; misunderstood the point of the game
                             mq.df.long$response.recode == "thats.the.only.way.it.can.come.out" | # Only + possibility modal
                             mq.df.long$response.recode == "theres.no.guarantee.so.i.guess.no" | # Adult; we interpret "guarantee" as a modal; understood modal; misunderstood the point of the game
                             mq.df.long$response.recode == "we.don't.know.maybe" | # Reading this as a correction
                             mq.df.long$response.recode == "yeah.it.doesnt.have.to.it.might.not.that.is"] <- "poss.modal" # Negated necessity is possibility

mq.df.long$response.recode[mq.df.long$response.recode == "i.dont.know" |
                             mq.df.long$response.recode == "i'm.not.sure" |
                             mq.df.long$response.recode == "im.not.so.sure" |
                             mq.df.long$response.recode == "no.or.not.sure" | # Correction
                             mq.df.long$response.recode == "not.sure.ill.say.no" | # Answer and then guess
                             mq.df.long$response.recode == "you.dont.know" |
                             mq.df.long$response.recode == "you.never.know" |
                             mq.df.long$response.recode == "i.don't.know"] <- "dont.know"

mq.df.long$response.recode[mq.df.long$response.recode == "i.hope.so " |
                             mq.df.long$response.recode == "i.hope.so" |
                             mq.df.long$response.recode == "it.will.but.it.will.also.come.out.over.here" | #OTHER
                             mq.df.long$response.recode == "first.time.it.will.come.out.here.second.time.it.will.come.out.here" |
                             mq.df.long$response.recode == "" |
                             mq.df.long$response.recode == "not.always" |
                             mq.df.long$response.recode == "pointed.to.uncertain.outside" |
                             mq.df.long$response.recode == "sometimes"] <- "other"

table(mq.df.long$response.recode)

mq.df.long$response.recode <- factor(mq.df.long$response.recode, levels = c("other", "dont.know", "negated.poss.modal", "poss.modal", "no", "yes"))

mq.df.long <- mq.df.long %>%
  mutate(discriminating = ifelse(input == "y" & output == "uncertain", 1, 0))

## Descriptive features of recoded responses ----

ggplot(data = mq.df.long, aes(x = age.group)) +
  geom_bar(aes(fill = response.recode), position = "fill") +
  facet_wrap(~ verb + input + output) +
  theme_bw() +
  theme(text = element_text(size=18), axis.ticks.y=element_blank(),
        panel.grid.major.x = element_blank(), panel.grid.minor.x = element_blank(), axis.ticks = element_blank(),
        panel.background = element_blank(), legend.title = element_blank(),
        strip.background = element_blank()) 

# Shortened version of the above figure with no "will" data

ggplot(data = mq.df.long[mq.df.long$verb != "will",], aes(x = age.group)) +
  geom_bar(aes(fill = response.recode), position = "fill") +
  facet_wrap(~ verb + input + output, nrow = 2) +
  theme_bw() +
  theme(text = element_text(size=18), axis.ticks.y=element_blank(),
        panel.grid.major.x = element_blank(), panel.grid.minor.x = element_blank(), axis.ticks = element_blank(),
        panel.background = element_blank(), legend.title = element_blank(),
        strip.background = element_blank()) 

# Number of each recoded response for each input, output, and age group
mq.response.recode.plot.table <- mq.df.long %>%
  group_by(age.group, verb, input, output) %>%
  summarize(num.yes = length(response.recode[response.recode == "yes"]),
            num.no = length(response.recode[response.recode == "no"]),
            num.poss.modal = length(response.recode[response.recode == "poss.modal"]),
            num.negated.poss.modal = length(response.recode[response.recode == "negated.poss.modal"]),
            num.dont.know = length(response.recode[response.recode == "dont.know"]),
            num.other = length(response.recode[response.recode == "other"]))

## Error analysis ----

table(mq.df.long$response.recode)


error.df <- mq.df.long %>%
  filter(verb != "will") %>%
  select(id, age.group, verb, input, output, response.recode) 

errors.coded <- error.df %>%
  group_by(id, age.group, verb) %>%
  summarize(strategy = ifelse(input == "straight" & output == "certain" & response.recode == "yes" &
                                input == "straight" & output == "uncertain" & (response.recode == "no" | response.recode == "negated.poss.modal") &
                                input == "y" & output == "certain" & (response.recode == "no" | response.recode == "negated.poss.modal") &
                                input == "y" & output == "uncertain" & (response.recode == "yes" | response.recode == "poss.modal"), "as.can",
                              ifelse(input == "straight" & output == "certain" & response.recode == "yes" &
                                       input == "straight" & output == "uncertain" & (response.recode == "no" | response.recode == "negated.poss.modal") &
                                       input == "y" & output == "certain" & (response.recode == "no" | response.recode == "negated.poss.modal") &
                                       input == "y" & output == "uncertain" & (response.recode == "no" | response.recode == "poss.modal"), "as.hafta", "doot")))

err.df <- mq.df.long %>%
  filter(verb != "will") %>%
  select(id, age.group, question, response.recode) %>%
  pivot_wider(names_from = question, values_from = response.recode) %>%
  mutate(can.strat = ifelse((can.straight.certain == "yes" | can.straight.certain == "poss.modal") & 
                              (can.straight.uncertain.in == "no" | can.straight.uncertain.in == "negated.poss.modal") & 
                              (can.straight.uncertain.out == "no" | can.straight.uncertain.out == "negated.poss.modal") &
                              (can.y.certain == "no" | can.y.certain == "negated.poss.modal") & 
                              (can.y.uncertain.in == "yes" | can.y.uncertain.in == "poss.modal") & 
                              (can.y.uncertain.out == "yes" | can.y.uncertain.out == "poss.modal"), "as.can",
                            ifelse(can.straight.certain == "yes" & can.straight.uncertain.in == "yes" & can.straight.uncertain.out == "yes" &
                                     can.y.certain == "yes" & can.y.uncertain.in == "yes" & can.y.uncertain.out == "yes", "yes.bias",
                                   ifelse((can.straight.certain == "yes" | can.straight.certain == "poss.modal") & 
                                            (can.straight.uncertain.in == "no" | can.straight.uncertain.in == "negated.poss.modal") & 
                                            (can.straight.uncertain.out == "no" | can.straight.uncertain.out == "negated.poss.modal") &
                                            (can.y.certain == "no" | can.y.certain == "negated.poss.modal") & 
                                            (can.y.uncertain.in == "no" | can.y.uncertain.in == "poss.modal") & 
                                            (can.y.uncertain.out == "no" | can.y.uncertain.out == "poss.modal"), "as.hafta",
                                          ifelse(can.straight.certain == "yes" & can.straight.uncertain.in == "no" & can.straight.uncertain.out == "no" &
                                                   can.y.certain == "no" & ((can.y.uncertain.in == "yes" & can.y.uncertain.out == "no") |
                                                                              (can.y.uncertain.in == "no" & can.y.uncertain.out == "yes")), "knows.which", "other")))),
         hafta.strat = ifelse((hafta.straight.certain == "yes" | hafta.straight.certain == "poss.modal") & 
                                (hafta.straight.uncertain.in == "no" | hafta.straight.uncertain.in == "negated.poss.modal") & 
                                (hafta.straight.uncertain.out == "no" | hafta.straight.uncertain.out == "negated.poss.modal") &
                                (hafta.y.certain == "no" | hafta.y.certain == "negated.poss.modal") & 
                                (hafta.y.uncertain.in == "yes" | hafta.y.uncertain.in == "poss.modal") & 
                                (hafta.y.uncertain.out == "yes" | hafta.y.uncertain.out == "poss.modal"), "as.can",
                              ifelse(hafta.straight.certain == "yes" & hafta.straight.uncertain.in == "yes" & hafta.straight.uncertain.out == "yes" &
                                       hafta.y.certain == "yes" & hafta.y.uncertain.in == "yes" & hafta.y.uncertain.out == "yes", "yes.bias",
                                     ifelse((hafta.straight.certain == "yes" | hafta.straight.certain == "poss.modal") & 
                                              (hafta.straight.uncertain.in == "no" | hafta.straight.uncertain.in == "negated.poss.modal") & 
                                              (hafta.straight.uncertain.out == "no" | hafta.straight.uncertain.out == "negated.poss.modal") &
                                              (hafta.y.certain == "no" | hafta.y.certain == "negated.poss.modal") & 
                                              (hafta.y.uncertain.in == "no" | hafta.y.uncertain.in == "poss.modal") & 
                                              (hafta.y.uncertain.out == "no" | hafta.y.uncertain.out == "poss.modal"), "as.hafta",
                                            ifelse(hafta.straight.certain == "other" | hafta.straight.uncertain.in == "other" | hafta.straight.uncertain.out == "other" |
                                                     hafta.y.certain == "other" | hafta.y.uncertain.in == "other" | hafta.y.uncertain.out == "other", "other",
                                                   ifelse(hafta.y.uncertain.out == "dont.know", "other",
                                                          ifelse(hafta.straight.uncertain.in == "poss.modal" | hafta.straight.uncertain.out == "poss.modal" |
                                                                   hafta.straight.certain == "dont.know", "other",
                                                                 ifelse(hafta.straight.certain == "yes" & hafta.straight.uncertain.in == "no" & hafta.straight.uncertain.out == "no" &
                                                                          hafta.y.certain == "no" & ((hafta.y.uncertain.in == "yes" & hafta.y.uncertain.out == "no") |
                                                                                                       (hafta.y.uncertain.in == "no" & hafta.y.uncertain.out == "yes")), "knows.which",
                                                                        ifelse(hafta.straight.certain == "no" & hafta.straight.uncertain.in == "no" & hafta.straight.uncertain.out == "no" &
                                                                                 hafta.y.certain == "no" & hafta.y.uncertain.in == "no" & hafta.y.uncertain.out == "no", "no.bias", "other")))))))))

err.df$overall.strat <- paste0(rep("CAN.", nrow(err.df)), err.df$can.strat, rep(".HAFTA.", nrow(err.df)), err.df$hafta.strat)

err.df.plot <- err.df %>%
  mutate(overall.strat = factor(overall.strat, levels = c("CAN.other.HAFTA.other",
                                                          "CAN.other.HAFTA.knows.which",
                                                          "CAN.other.HAFTA.as.can",
                                                          "CAN.as.can.HAFTA.other",
                                                          "CAN.yes.bias.HAFTA.yes.bias",
                                                          "CAN.yes.bias.HAFTA.as.can",
                                                          "CAN.as.can.HAFTA.yes.bias",
                                                          "CAN.as.hafta.HAFTA.no.bias",
                                                          "CAN.as.can.HAFTA.no.bias",
                                                          "CAN.knows.which.HAFTA.knows.which",
                                                          "CAN.as.can.HAFTA.knows.which",
                                                          "CAN.as.can.HAFTA.as.can",
                                                          "CAN.as.can.HAFTA.as.hafta"
  ),
  labels = c("Both other",
             "Can other, Have to knows which",
             "Can other, Have to as Can",
             "Can as can, Have to other",
             "Both yes bias",
             "Can yes bias, Have to as can",
             "Can as can, Have to yes bias",
             "Can as have to, Have to no bias",
             "Can as can, Have to no bias",
             "Can knows which, Have to knows which",
             "Can as can, Have to knows which",
             "Both as can",
             "Both correct")),
  overall.strat.simp = factor(ifelse(overall.strat == "Both correct", "Both correct",
                                     ifelse(overall.strat == "Both as can", "Both as can", "Other")), 
                              levels = c("Other", "Both as can", "Both correct")))

ggplot(err.df.plot[err.df.plot$overall.strat.simp != "Other",], aes(x = age.group)) +
  geom_bar(aes(fill = overall.strat.simp)) +
  geom_text(data = n.by.age, aes(x = age.group, y = 24.5, label = count.exp)) +
  scale_fill_manual(values=cbPalette[c(5,4)]) +
  theme_bw() +
  theme(text = element_text(size=18), axis.ticks.y=element_blank(),
        panel.grid.major.x = element_blank(), panel.grid.minor.x = element_blank(), panel.grid.minor.y = element_blank(),
        panel.background = element_blank(), legend.position = "",
        strip.background = element_blank()) +
  labs(x = "Age group", y = "Number of participants")

## Inferential analysis of responses on the linguistic task ----

# Do kids give "hafta"-appropriate responses more to "hafta" than they do to "can"?
# That is, is there any sign that of differentiation of the two verbs?
# We'll also do the converse: 
# Do kids give "Can-appropriate" responses more for "Can" than for "hafta"?

can.hafta.response.df <- mq.df.long %>%
  filter(input == "y" & output == "uncertain") %>%
  select(id, age.group, verb, response.recode) %>%
  mutate(hafta.appropriate.response = ifelse(response.recode == "no" | response.recode == "poss.modal", 1, 0),
         can.appropriate.response = ifelse(response.recode == "yes" | response.recode == "poss.modal", 1, 0))

hafta.response.mem <- glmer(hafta.appropriate.response ~  age.group * verb + (1 | id), data = can.hafta.response.df[can.hafta.response.df$verb != "will",],
                            family = binomial, control = glmerControl(optimizer = "optimx", calc.derivs = FALSE,
                                                                      optCtrl = list(method = "nlminb", starttests = FALSE, kkt = FALSE)))

hafta.response.emm <- emmeans(hafta.response.mem, specs = "age.group", by = c("age.group", "verb"), type = "response", infer = TRUE)
hafta.response.plot <- as.data.frame(hafta.response.emm)

ggplot(hafta.response.plot[hafta.response.plot$age.group != "adult",], aes(x = age.group, y = prob, group = verb, color = verb)) +
  geom_point(position = position_dodge(width = .35)) +
  geom_errorbar(aes(ymin = asymp.LCL, ymax = asymp.UCL), width = .05, position = position_dodge(width = .35)) +
  labs(title = "Probability of a hafta-appropriate response", x = "age group", y = "Model estimates and 95% Confidence Intervals") +
  theme_bw() +
  theme(text = element_text(size=18), axis.ticks.y=element_blank(),
        panel.grid.major.x = element_blank(), panel.grid.minor.x = element_blank(), axis.ticks = element_blank(),
        panel.background = element_blank(),
        strip.background = element_blank()) 

summary(hafta.response.mem)

pairs(emmeans(hafta.response.mem, specs = "verb", by = "age.group", type = "response", infer = TRUE))

can.response.mem <- glmer(can.appropriate.response ~  age.group * verb + (1 | id), data = can.hafta.response.df[can.hafta.response.df$verb != "will",],
                          family = binomial, control = glmerControl(optimizer = "optimx", calc.derivs = FALSE,
                                                                    optCtrl = list(method = "nlminb", starttests = FALSE, kkt = FALSE)))

can.response.emm <- emmeans(can.response.mem, specs = "age.group", by = c("age.group", "verb"), type = "response", infer = TRUE)
can.response.plot <- as.data.frame(can.response.emm)
can.response.plot$asymp.LCL <- ifelse(can.response.plot$asymp.LCL < .000001, NA, can.response.plot$asymp.LCL) # Eliminate full-panel error bars
can.response.plot$asymp.UCL <- ifelse(can.response.plot$asymp.LCL < .000001, NA, can.response.plot$asymp.UCL) # Eliminate full-panel error bars

ggplot(can.response.plot[can.response.plot$age.group != "adult",], aes(x = age.group, y = prob, group = verb, color = verb)) +
  geom_point(position = position_dodge(width = .35)) +
  geom_errorbar(aes(ymin = asymp.LCL, ymax = asymp.UCL), width = .05, position = position_dodge(width = .35)) +
  labs(title = "Probability of a can-appropriate response", x = "age group", y = "Model estimates and 95% Confidence Intervals") +
  theme_bw() +
  theme(text = element_text(size=18), axis.ticks.y=element_blank(),
        panel.grid.major.x = element_blank(), panel.grid.minor.x = element_blank(), axis.ticks = element_blank(),
        panel.background = element_blank(), 
        strip.background = element_blank()) 

summary(can.response.mem)

pairs(emmeans(can.response.mem, specs = "verb", by = "age.group", type = "response", infer = TRUE))

# Let's put 'can' and 'hafta' into one plot

hafta.response.plot.adj <- hafta.response.plot %>%
  filter(age.group != "adult") %>%
  mutate(response.type = "Appropriate for 'have to'",
         verb = ifelse(verb == "hafta", "have to", "can"))

can.hafta.response.plot <- can.response.plot %>%
  filter(age.group != "adult") %>%
  mutate(response.type = "Appropriate for 'can'",
         verb = ifelse(verb == "hafta", "have to", "can")) %>%
  rbind(., hafta.response.plot.adj) %>%
  mutate(response.type = factor(response.type, levels = c("Appropriate for 'have to'", "Appropriate for 'can'")))

ggplot(can.hafta.response.plot, aes(x = age.group, y = prob, group = verb, color = verb)) +
  facet_wrap(~response.type) +
  geom_point(position = position_dodge(width = .35)) +
  geom_errorbar(aes(ymin = asymp.LCL, ymax = asymp.UCL), width = .05, position = position_dodge(width = .35)) +
  labs(title = "Probability of an appropriate response", x = "age group", y = "Model estimates and 95% Confidence Intervals") +
  theme_bw() +
  theme(text = element_text(size=18), axis.ticks.y=element_blank(),
        panel.grid.major.x = element_blank(), panel.grid.minor.x = element_blank(), axis.ticks = element_blank(),
        panel.background = element_blank(), 
        strip.background = element_blank()) 

## Next test catching performance over trials ----

catching.mq.df <- read.csv("catching.mq.df.csv") %>%
  select(id, age.group, trial, response) %>%
  pivot_wider(names_from = "trial", values_from = "response") %>%
  mutate(can.correct.overall = ifelse(err.df$can.strat == "as.can", "correct", "incorrect"), # Use the error data to code whether they got "can" correct or not
         hafta.correct.overall = ifelse(err.df$hafta.strat == "as.hafta", "correct", "incorrect")) %>% # Ditto for whether they got "have to" correct or not
  pivot_longer(pretest1:posttest6, names_to = "trial", values_to = "response") %>%
  mutate(trial.type = ifelse(grepl("pre", trial), "pre",
                             ifelse(grepl("int", trial), "int", "pos")),
         trial.number = ifelse(trial.type == "pre", as.numeric(substr(trial, nchar(trial), nchar(trial))),
                                          ifelse(trial.type == "int", as.numeric(substr(trial, nchar(trial), nchar(trial))) + 4,
                                                 as.numeric(substr(trial, nchar(trial), nchar(trial))) + 10)),
         correct = ifelse(response == "certain", 1, 0))


catching.by.trial <- glmer(correct ~ trial.number * age.group + (1 | id), family = binomial, data = catching.mq.df,
                           control = glmerControl(optimizer = "optimx", calc.derivs = FALSE,
                                                  optCtrl = list(method = "nlminb", starttests = FALSE, kkt = FALSE)))

catching.by.trial.trends <- emtrends(catching.by.trial, specs = ~ 1 + age.group, var = "trial.number", infer = TRUE)
exp(as.data.frame(catching.by.trial.trends)[,c(2, 5, 6)])
catching.by.trial.trends.plot <- emmip(catching.by.trial, age.group ~ trial.number, cov.reduce = range, type = "response", plotit = FALSE)

ggplot(catching.by.trial.trends.plot, aes(x = xvar, y = yvar, color = age.group)) +
  geom_line() +
  ylim(0,1) +
  theme_bw() +
  theme(text = element_text(size=18), axis.ticks.y=element_blank(),
        panel.grid.minor.x = element_blank(), axis.ticks = element_blank(),
        panel.background = element_blank(), 
        strip.background = element_blank()) +
  scale_x_continuous(breaks = seq(0, 15, 5), labels=c(1, 6, 11, 16)) +
  labs(y = "Probability of placing \nwagon under target", x = "Trial number")

## Binomial tests of marble catching + Plot ----
binomial.analysis <- catching.mq.df %>%
  select(id, age.group, trial, response) %>%
  mutate(id = as.numeric(id)) %>%
  group_by(id, age.group) %>%
  summarise(num.correct = length(response[response == "certain" | response == "certan"]), 
            num.trials = length(response),
            binomial.p = binom.test(num.correct, num.trials, p = .5)$p.value,
            significant = ifelse(binomial.p < .05, "*", "ns"),
            significantly.better = factor(ifelse(binomial.p < .05 & num.correct > 8, "sig. better than chance", 
                                                 ifelse(binomial.p < .05 & num.correct < 8, "sig. worse than chance", "not sig.")), 
                                          levels = c("sig. better than chance", "not sig.", "sig. worse than chance")))

table(binomial.analysis$age.group)

binomial.analysis %>%
  group_by(age.group) %>%
  count(significantly.better)

ggplot(binomial.analysis) +
  geom_histogram(aes(x = num.correct, fill = significantly.better)) +
  facet_wrap(~age.group, nrow = 1) +
  theme_bw() +
  theme(text = element_text(size=22), panel.spacing = unit(2, "lines"),
        panel.grid.minor = element_blank(), panel.grid.major.x = element_blank(),
        legend.title = element_blank()) +
  labs(x = "Number of wagons under nonbranching slide", y = "Number of participants")

## Next test the relationship between performance on 'can' and catching ----

can.catching.by.test <- glmer(correct ~ age.group * trial.type * can.correct.overall + (1 | id), family = binomial, data = catching.mq.df[catching.mq.df$trial.type != "int",],
                              control = glmerControl(optimizer = "optimx", calc.derivs = FALSE,
                                                     optCtrl = list(method = "nlminb", starttests = FALSE, kkt = FALSE)))

summary(can.catching.by.test)
can.catching.by.test.emm <- emmeans(can.catching.by.test, specs = "trial.type", by = c("can.correct.overall", "age.group"), type = "response", infer = TRUE)
can.catching.by.test.emm.plot <- as.data.frame(can.catching.by.test.emm)
can.catching.by.test.emm.plot$age.group <- factor(can.catching.by.test.emm.plot$age.group, labels = c("4yos", "5yos", "6yos", "7yos", "adults"))
can.catching.by.test.emm.plot$can.correct.overall <- factor(can.catching.by.test.emm.plot$can.correct.overall, labels = c("can questions incorrect", "can questions correct"))
can.catching.by.test.emm.plot$trial.type <- factor(can.catching.by.test.emm.plot$trial.type, labels = c("before q's", "after q's"))
can.catching.by.test.emm.plot$asymp.LCL <- ifelse(can.catching.by.test.emm.plot$asymp.LCL < .00000001, NA, can.catching.by.test.emm.plot$asymp.LCL) # Remove full-panel error bars
can.catching.by.test.emm.plot$asymp.UCL <- ifelse(can.catching.by.test.emm.plot$asymp.LCL < .00000001, NA, can.catching.by.test.emm.plot$asymp.UCL) # Remove full-panel error bars

pairs(can.catching.by.test.emm) # Significant improvement among 4yos can-correct, but nowhere else

# When we collapse over age group and trial type, are kids who got can-questions wrong significantly better than .5?
# We need to refit the model, eliminating the 7yos and adults (due to NAs)
can.catching.simple <- glmer(correct ~ can.correct.overall + (1 | id), family = binomial, 
                             data = catching.mq.df[catching.mq.df$trial.type != "int" &
                                                     catching.mq.df$age.group != "7" &
                                                     catching.mq.df$age.group != "adult",],
                             control = glmerControl(optimizer = "optimx", calc.derivs = FALSE,
                                                    optCtrl = list(method = "nlminb", starttests = FALSE, kkt = FALSE)))
emmeans(can.catching.simple, specs = "can.correct.overall", type = "response", infer = TRUE)

# Let's add dotplots with individual proportions correct

prop.correct.dotplot.can <- catching.mq.df[catching.mq.df$trial.type != "int",] %>%
  group_by(id, age.group, can.correct.overall, trial.type) %>%
  summarize(prop.correct = mean(correct)) %>%
  mutate(trial.type = factor(trial.type, levels = c("pre", "pos"), labels = c("before q's", "after q's")),
         age.group = factor(age.group, levels = c("4", "5", "6", "7", "adult"),
                            labels = c("4yos", "5yos", "6yos", "7yos", "adults")),
         can.correct.overall = factor(can.correct.overall, levels = c("incorrect", "correct"),
                                      labels = c("can questions incorrect", "can questions correct")))

can.catching.plot <- ggplot(can.catching.by.test.emm.plot[can.catching.by.test.emm.plot$age.group != "adults",], aes(x = trial.type, y = prob)) +
  geom_point() +
  geom_errorbar(aes(ymin = asymp.LCL, ymax = asymp.UCL), width = .05) +
  geom_dotplot(data = prop.correct.dotplot.can[prop.correct.dotplot.can$age.group != "adults",], aes(y = prop.correct), binaxis = "y", position = position_nudge(x = -0.1),
               dotsize = 1, stackdir = "down", right = FALSE, alpha = .25) +
  facet_wrap(~ age.group + can.correct.overall, nrow = 5) +
  ylim(0,1) +
  theme_bw() +
  theme(text = element_text(size=18), axis.ticks.y=element_blank(),
        panel.grid.major.x = element_blank(), panel.grid.minor.x = element_blank(), axis.ticks = element_blank(),
        panel.background = element_blank(), 
        strip.background = element_blank()) +
  labs(y = "Probability of placing wagon under target branch", x = "Before vs. after modal questions")

can.catching.plot

## Next test the relationship between performance on 'hafta' and catching ----

hafta.catching.by.test <- glmer(correct ~ age.group * trial.type * hafta.correct.overall + (1 | id), family = binomial, data = catching.mq.df[catching.mq.df$trial.type != "int" & catching.mq.df$id != 509,],
                                control = glmerControl(optimizer = "optimx", calc.derivs = FALSE, # Participant 509 needs to be categorized for "hafta.correct.overall"
                                                       optCtrl = list(method = "nlminb", starttests = FALSE, kkt = FALSE)))

summary(hafta.catching.by.test)

pairs(emmeans(hafta.catching.by.test, specs = "hafta.correct.overall", by = "age.group"), type = "response") # No significant results for the effect of getting the questions right at any age group, collapsing over trial type
pairs(emmeans(hafta.catching.by.test, specs = "age.group")) # No significant differences between any age groups, collapsing over trial type and whether they got the hafta questions correct
emmeans(hafta.catching.by.test, specs = "age.group", by = "hafta.correct.overall", type = "response", infer = TRUE)
pairs(emmeans(hafta.catching.by.test, specs = "trial.type", by = c("hafta.correct.overall", "age.group"), type = "response", infer = TRUE))

hafta.catching.by.test.emm <- emmeans(hafta.catching.by.test, specs = "trial.type", by = c("hafta.correct.overall", "age.group"), type = "response", infer = TRUE)
hafta.catching.by.test.emm.plot <- as.data.frame(hafta.catching.by.test.emm)
hafta.catching.by.test.emm.plot$age.group <- factor(hafta.catching.by.test.emm.plot$age.group, labels = c("4yos", "5yos", "6yos", "7yos", "adults"))
hafta.catching.by.test.emm.plot$hafta.correct.overall <- factor(hafta.catching.by.test.emm.plot$hafta.correct.overall, labels = c("'have to' questions incorrect", "'have to' questions correct"))
hafta.catching.by.test.emm.plot$trial.type <- factor(hafta.catching.by.test.emm.plot$trial.type, labels = c("before q's", "after q's"))

hafta.catching.by.test.emm.plot$asymp.LCL <- ifelse(hafta.catching.by.test.emm.plot$asymp.LCL < .00000001, NA, hafta.catching.by.test.emm.plot$asymp.LCL) # Remove full-panel error bars
hafta.catching.by.test.emm.plot$asymp.UCL <- ifelse(hafta.catching.by.test.emm.plot$asymp.LCL < .00000001, NA, hafta.catching.by.test.emm.plot$asymp.UCL) # Remove full-panel error bars

# Let's add dotplots with individual proportions correct

prop.correct.dotplot.hafta <- catching.mq.df[catching.mq.df$trial.type != "int",] %>%
  group_by(id, age.group, hafta.correct.overall, trial.type) %>%
  summarize(prop.correct = mean(correct)) %>%
  mutate(trial.type = factor(trial.type, levels = c("pre", "pos"), labels = c("before q's", "after q's")),
         age.group = factor(age.group, levels = c("4", "5", "6", "7", "adult"),
                            labels = c("4yos", "5yos", "6yos", "7yos", "adults")),
         hafta.correct.overall = factor(hafta.correct.overall, levels = c("incorrect", "correct"),
                                        labels = c("'have to' questions incorrect", "'have to' questions correct"))) %>%
  filter(!is.na(hafta.correct.overall) & age.group != "adults")

hafta.catching.by.test.emm.plot <- hafta.catching.by.test.emm.plot %>%
  filter(age.group != "adults") 

hafta.catching.plot <- ggplot(hafta.catching.by.test.emm.plot, aes(x = trial.type, y = prob)) +
  geom_point() +
  geom_errorbar(aes(ymin = asymp.LCL, ymax = asymp.UCL), width = .05) +
  geom_dotplot(data = prop.correct.dotplot.hafta, aes(y = prop.correct), binaxis = "y", position = position_nudge(x = -0.1),
               dotsize = 1, stackdir = "down", right = FALSE, alpha = .25) +
  facet_wrap(~ age.group + hafta.correct.overall, nrow = 5) +
  ylim(0,1) +
  theme_bw() +
  theme(text = element_text(size=18), axis.ticks.y=element_blank(),
        panel.grid.major.x = element_blank(), panel.grid.minor.x = element_blank(), axis.ticks = element_blank(),
        panel.background = element_blank(), 
        strip.background = element_blank()) +
  labs(y = "Probability of placing wagon under target branch", x = "Before vs. after modal questions")
hafta.catching.plot

# I'm going to make these 2 figures into 1

ggarrange(plots = list(can.catching.plot, hafta.catching.plot), nrow = 1)
